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1. Introduction 

Two regimes are often distinguished when solving 
structural dynamics and vibroacoustics problems in the 
frequency domain. Broadly speaking, 

(i) the low-frequency range is a region with low 
modal density and low response sensitivity, 
whereas 

(ii) the high-frequency range is a region with high 
modal and statistical overlaps. 

The reasons for this distinction are multiple, with four 
salient characteristics: all details of the behaviour of the 
systems are not required or requested from a practical 
engineering point of view; the cost of modelling and 
solving a system in all its complexity may be excessive; 
some parts of the acoustic or vibratory response may 
be so sensitive to the system parameters that their 
adequacy to represent the system behaviours precisely is 
at best tenuous at some frequencies; and, finally, a limited 
set of simplified features or asymptotic behaviours 
may actually dominate the response in some frequency 
ranges. 

© 2014 The Authors. Published by the Royal Society under the terms of the 
Creative Commons Attribution License http://creativecommons.org/licenses/ 
by/3.0/, which permits unrestricted use, provided the original author and 
source are credited. 



In the low-frequency range, a relatively small number of well-defined modes are usually 
sufficient to describe the response of the system deterministically Furthermore, thanks to 
long wavelengths of its mode shapes and propagating waves, element-based models of the 
system can be kept coarse, and therefore at a reasonable size, while still providing accurate 
representation of its behaviour. Usually, the system response is therefore essentially treated in 
a deterministic manner at low frequency Relatively unexpensive physical models give accurate 
representation of the system behaviour and a small number of modes or waves can be evaluated 
and used. 

The same approach can usually not be used at high frequency. Dimensionally huge models 
are necessary owing to small wavelengths, the number of modes can be extremely large, and 
treating responses deterministically does not really make sense owing to their sensitivity. In this 
frequency range, statistical or average methods are usually preferred. An additional motivating 
factor for this choice, besides the facts that a full response could be very expensive and practically 
random in the presence of small system disturbances or propagating numerical errors, is that the 
statistical properties of the response may greatly simplify. 

This fact has been most notably used in statistical energy analysis (SEA) [1,2], a now common 
branch of structural dynamics and vibro-acoustics whose developments were sparked by Lyon 
and Maidanik's paper 50 years ago [3]. Under certain conditions [4], the response of structural 
or acoustic components can indeed be characterized by their energy density in a frequency band 
and the exchange of energy between these components is proportional to the difference between 
their energy densities. This simplification has been used successfully to evaluate the vibration 
and acoustics behaviour of complex systems such as ships whose response would otherwise have 
been mostly untractable. A library of SEA parameters, such as coupling loss factors, has been 
developed for a range of material components and geometries in an analogy to the library of 
deterministic components developed at low frequency. Progress in the understanding of SEA and 
trying to extend its applicability has been made [5-8] but without generally allowing it to overlap 
fully with the lower frequency methods. 

This problematic mid-frequency gap in understanding and modelling between low- and high- 
frequency regions has been long known and an area of concern. While there is an increasingly 
strong industrial and theoretical call for techniques that can be applied in the whole frequency 
range, and while both the physical or modal lower frequency techniques and the higher frequency 
SEA methods are pushed as involuntary candidate to fill this gap, neither approach appears to be 
extendable enough to cover the whole mid-frequency region. An additional difficulty is that the 
boundaries of applicability of low- or high-frequency methods are somewhat blurred and it may 
thus be hard to assess whether those are adequate in particular engineering situations. 

There have been extensive efforts to describe and extend the range of applicability of the 
respective methods such as, starting from the higher frequencies point of view, in [9-15]. 
Similar efforts to encroach on the mid-frequency range have been pursued, starting from the 
lower frequencies point of view. For example, wave methods that allow to obtain deterministic 
responses over large physical domains using only a relatively small number of waves have been 
further developed. Since the wave characteristics can be evaluated from a detailed model of only 
a small portion of the domain, it is possible to part with several limitations of low-frequency 
methods. Alternative statistical simplifications at high frequencies have also been considered such 
as the representation of the system dynamic matrices by random matrices. A description of some 
recent developments can be found in [16-18]. 

Significant advances have also particularly been made in the combination of deterministic 
and statistical features of system characteristics. Besides the techniques to include parametric 
uncertainty in the low-frequency range, most notable, are relatively recent hybrid techniques such 
as those proposed in [19-24] that specifically deal with the mid-frequency range. The argument 
is that in the same structure or system, different components may see a given frequency as being 
either in the low- or in the high-frequency ranges. The latter 'fuzzy', 'complex' or 'reverberant' 
components are coupled to the former deterministic 'master' deterministic components, so 



that displacement variables are interfaced with energy variables and vice-versa. This allows 
application to specific practical engineering problems that would only be treatable with great 
difficulty otherwise, for example in the case of car doors that can be seen as made of relatively 
low-frequency frames coupled to very thin, high-frequency plates. 

It is the view of the author that there remains a need for general points of view in which 
the various methods can be assessed, compared and extended. The work presented here intends 
to offer such a point of view that covers the whole frequency range. The proposed approach 
is to consider averaging or filtering of the frequency responses with tunable averaging width. On 
the one hand, a small averaging width leads to dealing with the system deterministically, with the 
frequency average being asymptotically equal to the deterministic response and zero variance. On 
the other hand, a larger averaging width, leads to a statistical and energy treatment of the system 
behaviour. Since this averaging width can be freely chosen and tuned at different values along 
the frequency range, it is possible to use the same analysis with a smooth transition from the low- 
to the mid- and high-frequency ranges. 

This paper is organized in the following way. In §2, the proposed averaging process is 
formally described in general terms. It is then shown in §3 that the frequency average, covariance 
and averaged square of the response are more than theoretical concepts and that they can be 
evaluated in practical situations. The presentation focuses especially on the case of Gaussian 
averages and other possible types of averaging functions are mentioned. Illustrations are then 
presented and discussed in §4. Particular focus then turns to the time analysis of the averaged 
response in §5, and the preservation of the system energy information in §6. Finally, the 
efficient evaluation of the average response through matrix operations, without modal analysis, 
is mentioned in §7 before conclusions are drawn. 

The modal expressions of the averages are very general. For any system that can be expressed 
in the form of equation (3.2), the frequency average of its response has the form of (3.8), where 
the function S(-) takes the particular values of (3.11)-(3.13) in the case of a Gaussian average. 
Similarly the simple expressions of the variance, covariance, and the expectation of the squares 
or products of the responses can be found in (3.16) and (3.28), with the values of the function Q(-) 
being discussed in §3c. 



2. Proposed framework 



The approach proposed in this paper is applicable to linear systems such as 



A(co)x(co) = f. 



(2.1) 



where the relationships between the input and output vectors, f and x are described by a dynamic 
stiffness matrix A. This matrix and the response depend on a circular frequency parameter, co = 
Inf and the interest is in the response x{o)) = A{(jL>)~^f{oj), in a frequency range co e [(Wmin/<^max]- 
The full response at every possible frequency within the range is however usually not required. 
An engineer may indeed only be interested in one or more particular transfer functions, such 
as g(co) = c^x(co), for an output vector c, or some combination thereof. Or, most important in the 
context of this article, one may only have broad interest in the frequency range: one may only be 
interested in a few particular modes, in the value of the response in some discrete regions, or in 
the average of the response magnitude or energy. The analysis in the frequency domain can be 
transposed in the time (t) domain, through the Fourier transforms 



x((jo) e^^^ dco and inversely, x(co) = 



x(f)e-'^^ dt. 



(2.2) 



(a) Frequency average, averaged square and variance as computational goals 



Instead of considering either the deterministic solution or the energy statistics of the system 
described by equation (2.1) as is typically done at low- or high-frequency ranges, it is proposed 
here to directly consider the statistics or average of the response x{a)) with regards to a variation in 
the frequency The problem or computational goal is thus as follows: for a frequency-dependent 
level of uncertainty or averaging width, a{a)), and a frequency probability density function or 
frequency averaging function pai^oj), evaluate the average and variance, that is 

x(co) = E[x{(D + Aoj)] and var[x((X))] = E[abs(x(a; + Aoj) - x(co)f] (2.3) 

in the frequency range co, coa, o)b ^ [<^min/ <^max]- The expectation or average E[.] is on the random 
variable Aco so that, for any function u, 



E[u{(jL> + Aoj)] ■- 



u(a) + Aa))pa{^o)) dAco, (2.4) 



where Vp denotes the support of the measure Pa{-)f generally [—00,00] for pa{.) defined on 
the real axis. Similarly, the covariance is considered for possibly different input vectors and 
different frequencies, both distinguished by subindices .a and .5 so that A(a;A)xA(<^A) = 
and A((Wb)xb((Wb) = fs- It is defined by 

cov[xa((Wa), ^b{o^b)] = E[(xA(a;A + Ao;) - xa('^a))(xb(<^b + Ao;) - x^{ojb))^1 (2.5) 

where the superscript indicates the complex transpose conjugate. The variance of a vector is 
thus its covariance with itself. In general, when evaluating a covariance at different frequencies, 
the averaging width, a{(jOA, <^b)/ rnay depend on the two frequencies. 

The frequency average and covariance can be seen as the first and the second frequency 
moments of the response. However, since the responses are complex functions, some second 
moment information on the phase would be lost if only the covariance was considered. In order 
to remediate this, one also targets as computational objective the frequency average of the square 
of the response, that is 

squ[x(a;)] = E[x(a; + Aoj)x{co + Acofl (2.6) 

where the superscript J denotes the transposed vector. An additional objective, similar to the 
covariance, is then also the expected value of the product of responses at different frequencies or 
even due to different forces. It is defined as 

csq[xA(<^A)/ xb(<^b)] = E[xa(<^a + Aco)xb{cob + Acoj^l (2.7) 

where the averaging width, a(coA, <^b)/ rnay again depend on two parameters. 

Information on the phase is then retained thanks to the transpose — rather than the complex 
conjugate transpose — of XB(a;B) being also considered. The frequency /zrsf moments or averages of 
the real and imaginary parts of the response are indeed 

E[m{x{co + Aco))\ = m{x{a))) and E[^{x{a) + Aco))] = S(x(a;)), (2.8) 
and their second moments for consistent averaging widths can be derived as 

E[m{xA{coA + Aa;))^(xB(a;B + Acoj^)] = ^^(S + P), (2.9) 

E[?il{xA(oJA + Aa))p{xB{ojB + Aojj^)] = ^^(S - P), (2.10) 

E[S(xA(a;A + Aoj)Mxb{ojb + Aojj^)] = \^{S + P) (2.11) 

and E[S(xA(a;A + Aoj)P{xb{ojb + Aojf)] = - ^^(S - P), (2.12) 

where S = csq[xA((^A)/XB(<^B)] and P = E[xa((^a + Aa;)xB((X)B + Ao;)^] = cov[xa((^a)/ ^b{o^b)]-^ 
xa(<^a)xb(<^b)^. 



(b) Link to deterministic response and statistical energy analysis 



The main advantages of the proposed solution objective, in regards to providing a consistent 
framework from low- to high-frequency ranges, are (i) that it can be applied in the whole 
frequency range and (ii) that it covers usual deterministic and statistical energy approaches as 
particular or asymptotic cases. The deterministic, low-frequency approach indeed corresponds 
to the asymptotic case where the averaging width tends to zero 



while the power frequency average, can be expressed in terms of the average and variance using 



The statistical simplifications encountered at high-frequency, say in the SEA theory, can thus be 
examined in this context. Since the averaging width a{(jo) can depend on frequency, a smooth 
transition can be obtained through and between regions usually categorized as low-, mid- and 
high-frequency ranges. Such categorization is not necessary when both the first and the second 
frequency moments of the response are considered, although it may be very useful to make use 
of accurate approximations in the respective frequency ranges. 



In order to take full advantage of the proposed framework and to be able to use it in practical 
situations, frequency average, variance, covariance, averaged square and averaged product of 
responses should ideally have analytical expressions that can be evaluated accurately, robustly 
and cheaply. 

Such exact analytical expressions of average, variance, covariance and therefore of power 
average have been presented in [25] and further discussed in [26-28]. For a real Gaussian function 
p^(.) = Pa {.), both average and variance of any transfer function can be expressed in terms of the 
Faddeeva function w{z) = (1 + (H/s/tt) Jq e^ dt) [29], with different expressions depending 
on the imaginary parts of the system eigenvalues. In the case of purely real eigenvalues, the 
corresponding modal components of the averages have to be understood in a principal value 
sense as noted by Gautschi [30, p. 188]. 

It is worth stressing the important subtlety that, as a function of its eigenvalue, each modal 
component corresponds to different analytic functions in the two half-spaces of positive and 
negative part of this eigenvalue. Since these functions are neither analytic continuation of each 
other nor equal to the principal value function on the real axis, it is essential to understand what 
is meant by an undamped modal component. If one means a modal component that has negligible 
damping, then one should use the expression corresponding to the limit of vanishingly small 
positive damping. On the other hand, in the case of an actual undamped modal component, there 
is no other alternative, without further information, than to use the principal value expression. 
Physically, if the average is evaluated through Monte Carlo sampling, in this case, there will 
always be samples drawn close enough to the resonance for the estimated average to not converge. 
However, if these rogue samples that are within a radius, say p, of the resonance are excluded 
when estimating the average, there will be convergence and the converged value will tend to the 
predicted principal value for vanishing values of p. 

It was further demonstrated that all expressions could be evaluated efficiently and robustly 
in the whole frequency range. The theory was also extended to the case of complex normal 
variables [27]. In this paper, particular focus is put on the real Gaussian distribution, denoted 
by a superscript for real argument Aco, zero mean and standard deviation a{a)) or a{a)Af^B)f 
so that one has 



x{co) = lim x{co) and lim var[x(a;)] = 0, 



(2.13) 



E[x{co + Aco)x(co + Acof^] = vaY[x(co)] + x(co)x(co)^. 



(2.14) 



3. Analytical frequency average, variance and averaged square 




(3.1) 



Other options of averaging functions that readily fit in the proposed framework are mentioned in 
§3e. The presentation uses a modal point of view in order to retrieve the exact scalar integrals 
that were presented in [25,27]. The resulting average expressions may however be expressed 
in matrix form without requiring any explicit evaluation of eigenvalues or eigenvectors. This 
has been notably demonstrated in [31] where it was further shown that Krylov methods can be 
used to evaluate the average in a frequency range extremely efficiently. The existence of such 
efficient matrix evaluations is essential for practical application of the theory on dimensionally 
large discretized systems. 



(a) Frequency modal expression of the response 

For facility of presentation, it is assumed that equation (2.1) can be expressed in an equivalent 
form where the dynamic stiffness matrix is a linear function of the frequency parameter, that is 



(A«-c.A®)x(l)H = f« 



(3.2) 



with constant Aq ^ and A^^^ . This is always the case if, say, the original dynamic stiffness matrix, 
A{co), is a polynomial function of oj. For example, in the case of a damped structural system, 
this matrix may be a quadratic function of oj as A{(jl>) = (K + icoC — o/M) with constant stiffness, 
damping and mass matrices, K, C and M. The matrices, output and input vectors of the equivalent 
form (3.2) could then be 



K 


c" 




"o 


-iM" 




x{a)) 


and f(^^ = 


"f" 












0 


I 


il 


0 


ia)x{a)) 




0 



(3.3) 



The eigenvalues ojj and left- and right-eigenvectors, f^^^ and 0^^^ ^ 0, of the system then satisfy 



[ffViAf-cojAf) = 0 and (A« - a;yA«)0f = 0. 



(3.4) 



The matrix A^^^ is assumed regular and if the system further has a full set of eigenvectors, 
<t>j\i = 1, . . . ,N, a modal decomposition can be obtained by using the corresponding basis, 
0(1) = {<t>^ ... <t>^\. Binormalizing the left- and right-eigenvectors with respect to the A^^^ matrix, 
so that «J^(1)^A^/^0(1) = I, ^^(1)^a[J^0(1) = diag(a)y), and expressing x^^\cd) = 0^^^y(co), equation (3.2) 
can indeed be rewritten with diagonal matrix, as (diag(a;y) — col)y{co) = which results in 

the following modal expression of the response: 



Pl0f([^f]^f) 

{ojj - 0)) 



for eigenvalues and modes that satisfy 



[ff]TA(a;y) = 0, A(a;y)0j" = 0. 



(3.5) 



(3.6) 



The matrix Px denotes the projection matrix that allows to extract the response vector x{(jo) from 
the longer vector x^^\a)), i.e. x{oj) = F^x^^\oj). In the case of the damped example of equation 
(3.3), it is equal to Px = [I 0]. While the eigenvectors, 0y = Px0y^^ are therefore the first half 

of the linearized eigenvectors 0^^^ in the present example, the same expression can be found — 
and the theory derived below can be used — in more general cases. Note also that if the force 
vector was frequency dependent and could be expressed in rational form of co, then the present 
approach could still be used with expansions similar to (3.5) and an alternative meaning 
for the poles. 



(b) Frequency average 

In modal form, the frequency average of the response is thus the function 



j=l,...,N I 



(ojj — CO — u) 



Vaiu) du 



= Yl {<l>jifj'^f)S{o^,a,coj)}. 

j=l,...,N 



(3.7) 
(3.8) 



Ideally, the scalar integrals S{co,a,a)j) = jj^^{l/{a)j — co — u))pa{u) du or their combination 
necessitate computable exact expressions or accurate approximations which are discussed next. 
As discussed in [28], the effect of the averaging operation on the individual eigenvalues can be 
combined into the form of an (exact) averaging matrix H{(jl>), so that x(co) = H(a;)f, where 



i=i,...,N 



(3.9) 



The average of other transfer functions such as i(jox{co) that corresponds to the system's velocity 
can be obtained directly from the application of the theory to the system written in the form (3.4), 
with matrices and vectors defined as in equation (3.3). The full averaging matrix is then simply 



j=l,...,N 



(3.10) 



with the same functions S{co,a,coj) and the vector of averaged response (and, possibly, of its 

averaged transfer functions of time derivatives) is y!^\co) = H^^\co)f^^\ 

Real Gaussian case. In the case of a real Gaussian averaging function, averaging corresponds to 
applying a Gaussian filter to the response or evaluating its Weierstrass transform (see for example 
[32,33] and the references therein). 

Since oj is real, the exact expressions of the integrals are [25,29,34] 



S^^\co,a,coj) - 



— 1 7T 

' a{co)i 2 
i /tT 



a{co) 
.PV ~ 



/ 2 
i [jf 



a{a))y 2 



\V2a{aj)J ^ 

M ifS(a;v)<0 

{a)j - cof 



^a{co) 



- exp 



la{a))- 



(3.11) 
(3.12) 
(3.13) 



where the complex conjugate, denoted as .*, of both the Faddeeva function and the eigenvalue in 
the second equation are taken and where the third equation has to be understood in a principal 
value sense. For example, if all the eigenvalues have a positive imaginary part, expression (3.8) of 
the exact average becomes 



a(co) V 2 



y r0.(f .Tf)n; f 



(3.14) 



Such expression can also be expressed and evaluated directly in terms of transfer functions, the 
average response matrix {—i/a{co))-s/7T/2 |Xl;=i n[^;(^;^)^(('^; ~ <^)/(V2fl((X»)))]| or the system 
matrices, without requiring a modal decomposition, as discussed in [31]. 



(c) Frequency variance and covariance 

The frequency variance or covariance of the response is also available, from equation (2.14) and 
the following frequency cross-power average 



k=l,...,Ni=l,...,N ^ 



{coj -coA-u) {col -^B-u) 



Va{u) du 



= E {<l>j<l>f{fj^iA){fk^hTQ{coA,oJB,a,ojj,cok)}. 

i,k=l,...,N 



(3.15) 
(3.16) 



The integrals Q{o)a, <^B/ ^/ <^;/ <^k) = ix) ^^j ~ ^A — /{col ~ ~ ^))pM this 
expression have one of two forms depending if coj — coa and {co]^ — oj^T = col — co^ are equal 
or not [25]: 



if coj — COA ^0^1 — oJBf then, by considering partial fraction decomposition, 

S{(L>Af(l{o^Af(^B)fOJj) - S{(L>Bf(l{o^Af(^B)f(^l) 



Q{o)A,m,ci,o)j,o)k) = 
otherwise, if coj — coa = ^1 — <^B/ 



col - COB- + '^A 



Q{0JA, COB, a, 0)j, Ok = 0)^ - 0)A + cjob) = 



. Dp {o)j — COA ~ uf' 
and, in the real case, one has by integration by parts, 

1 dpa{u) 



Q{a)A, COB, ci, coj, a)k = (oJ - (OA-\-coB) = - 



{coj — COA ~ u) du 



du. 



(3.17) 



pa{u)du (3.18) 



(3.19) 



The frequency average power, cross-power, variance and covariance are thus directly available 
from expressions (3.11)-(3.13) in the first case where coj — coa col ~ ^B- 

Real Gaussian case. Considering a real Gaussian and in the first case where coj — coa col ~ ^B, 
one finds from equations (3.17) and (3.11)-(3.13) that attention must be paid to the signs of both 
coj and col. example, if the system is damped and all eigenvalues coj have positive imaginary 
part, their complex conjugate, col negative imaginary part, so that 



E^^\xa{coa + Aa;)xB(a;B + Ao;)"] 



a{o)A,coBn 2 



p f 

/ — X < 


E E " 




=1,...,]V;=1,...,N 



<t>^<t>f{f^iA){fk'iBT 



w{{ojj - 0JA)/{V2a{0JA.0JB))) + w*((a;fc - ojB)/{V2a{ojA,a)B))) 
col- OJB- + 



(3.20) 



Using equations (2.14) and (3.14), the covariance is then found to be 



—1 



J2 E <t>j<t>f{fi'iA){fk^hT 

k=l,...,Nj=l,...,N 



a{a)A,o)v>n 2 

^w{{o)j - a)A)/{^a{a)Arm))) + w*((a;fc - a)i^)/{^a{a)A,o)B))) 



OJj. - OJB - OJj + OJA 



a{ojA,ojB)y 2 VV2^ 



(3.21) 



Note that in the same example case of coj with positive imaginary part, if the further condition 
coj — coA = oJk — OJB holds, then the double pole integrals equal 



a^(coj — CO a) V 2 



71 / 9^(a;v — CO a) ^(<^7 — <^a)\ 
V2fl V2fl / 



where K{., .) denotes the Voigt function [35], defined as the real part of the Faddeeva function 
K{x,y) = ?ii[w(x + iy)] = (y/n) J!^^ exp(-fi)/{(x - tf + y^) d^ for real x and y. 

Alternatively, when (Wy — coa = o)1 — <^B/ the integral Q^^\coA,oj^,ci,ojj,co]^ = a;* — o^a + (^b) of 

equation (3.17) has a simple expression since dp^^\u)/du = —(u/a(coA,ci>^)^)p^^\u). This leads 
indeed to the following: 



Q^§^ (co, a, coj, co]^ = co* -coa + cob) = -^ 



{cOj — COA — u 

OJj - OA 



1 COi — COA I \ 



(3.23) 



p-^\u)Au (3.24) 



(3.25) 



and expressions (3.11)-(3.13) can again be used, depending on the sign of ^{coj). 



(d) Frequency averaged square and averaged product of responses 

The frequency square or averaged product of the response vectors are similarly available from 



csq[xA(a;A), Mm)] = XI 1 ^i^k(^i^^A){fk^h) 

k=l,...,N i=l,...,N i 



(coj — COA ~ u) (CO]^ — co^ — u) 



Va{u) du 



i,k=l,...,N 



(3.26) 
(3.27) 



The integrals R{coA,(OB,ci,(jOj,co]^) = jjy^{l/{coj — coa — u)){l/{co]^ — co^ — u))pa{u) du are thus simply 
R{coA,oj^,ci,ojj,cok) = Q{coAfO)Bf(^f(^jf(^^) and, following the explanation of the previous section, 
therefore also have one of two forms depending if coj — coa and co]^ — co^ are equal or not. 

Real Gaussian case. Considering again the real Gaussian case and the particular example where 
the system is damped and all eigenvalues coj have positive imaginary part, the first case to 



consider is coj — oja ^co^ — oj^ when 



ITT 

X 



fc=l,...,N;=l,...,N '- 

<^fc — <^B — + <^A 



. (3.28) 



In the general real Gaussian case, one has, when coj — coa = co]^ — co^ and for real coa and co^, that 
R^^\co, a, coj, cok = ojj — coa + <^b) = Q^^\oj, a, coj, co^ = to* — o^a + cob), whose expression is provided 
in equation (3.25). 



(e) Other distributions 

The general expressions presented in this paper can be directly applied to a wide variety of 
averaging, filter or distribution functions. Many choices of such functions can be found in various 
fields such as in the theory of filters [36] or in statistics [37]. One of the most notable among those 
might be the uniform distribution since it has been extensively used in acoustics and vibration 
to model the energetic response in frequency bands or octaves. Some explicit expressions of 
the necessary integrals can be adapted from the expressions of the frequency average of the 
energy given in [38, appendix] for the particular case of a proportionally damped system, i.e. 
when the complex poles come in particular pair combinations. Regarding the first frequency 
moment, i.e. the frequency average, an advantage of distributions or averaging functions that 
are in rational form, such as the Cauchy, Fisher's f. Student's t distributions, or the Bessel, 
Butterworth, Chebyshev filters, is that they may offer the computational advantage that the 
average can easily be expressed in matrix form, as a function of system solutions only. This 
can for example be used, together with the fact that the averaged input energy is proportional 
to the average transfer function, to approximate the averaged energy in a frequency band 
as in [39]. 

A more general approach for the average, proposed in [31] and inspired by Weideman's work 
on scalar rational approximations of the complex error function [40], allows working with rational 
matrix operations and by-passing any modal analysis, even for general non-rational distributions. 
Further using properties of Krylov subspaces, the main overall cost to evaluate the averages over 
a range a frequencies is only a few solutions, as demonstrated in [31] and further discussed in 
§8. A significant advantage worth noting is that the average at a particular frequency can be 
obtained with solutions at a single shifted frequency which corresponds to working with a system 
with additional damping. Any available factorization of the system dynamic stiffness matrix can 
thus be re-used. 

The choice of which frequency averaging function to use is dependent on a user's preference 
and on the objectives of the application of interest. The influence of a particular choice can 
further be quantified, notably by examining the resulting frequency moments in the time 
domain. This is discussed in %%Sa, 6 and 7. It is seen that for a constant averaging width, 
the averaging process results in the scaling of the exact impulse responses by the inverse 
Fourier transform of the averaging function. In statistics, this Fourier transform of a distribution 
function is called its characteristic function and is therefore readily available for standard 
distribution functions. 

Complex averaging windows are also generally supported. For example, the integrals 
necessary for a complex Gaussian were presented in [27]. Some care must be taken since the 
covariance may exhibit local singularities. 
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Figure 1. Illustration of the complex spring-mass system benchmark. 
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4. Illustration of the evaluation of the frequency averages 

The use of the general notions of frequency average and frequency variance can be geared 
towards specific applications, depending on the computational goals pursued. Various uses are 
illustrated in the next sections. In the current section, a benchmark problem that exhibits a mix of 
low- to higher frequency characteristics is first introduced. The average, variance and averaged 
square are then evaluated and discussed, in comparison to the system full solution. 



(a) Description of the benchmark problem 

The benchmark system considered here is made of a large oscillator connected to a set of 
substructures, as illustrated in figure 1. Similar prototypes of complicated systems made of a 
large mass with many attached spring-mass systems have recurrently appeared in the literature 
of the last 20 years. One of the first occurrences was, for example, the undamped system that 
Weaver [41] used to illustrate the (ensemble) averaging approach he proposed as an alternative 
and extension to previous work on fuzzy and complex structures as by Soize [20] and Pierce 
et al. [21]. Although apparently simple, such complex systems are representative of engineering 
systems that are conceptually difficult to model because they exhibit a combination of low-, mid- 
and high-frequency characteristics. Their behaviour having actually proved to be far from simple 
and rather particularly interesting, they have been used in a form or another either to illustrate 
various theories or to study them in their own right. A salient example of the latter is the analysis 
in [19,42-44] and other publications of how the undamped additional spring-masses provide 
apparent equivalent damping to the larger mass. 

Similar line of thought and motivation as that of others has been followed to select the current 
benchmark. Its parameters are, in particular, inspired from the first example of [19, section 7.1]. 
Two aspects however somewhat distinguish the current work from previous work on this type of 
problem. First, somewhat unusual is that the presence of both significant or smaller values of 
damping are supported. Second, more generally and arguably more importantly, the analysis 
is exact and general, so that neither approximations, such as for example an ad hoc smoothing 
of the modes, nor conditions on the distribution of the modes are necessary. While the chosen 
benchmark has damping, the proposed exact theory can also be applied without approximation 
to either undamped or very slightly damped systems (doing so, the remarks made in §3 about 
principal values must be kept in mind). It is finally worth reminding here that the averaging in the 
current work is on a property (the frequency) of an actual system, rather than over an ensemble 
of similar structures. 

The details of the model are now described and the frequency average and variance are 
presented for varying values of the tuning parameter, a{co). Units notations are essentially omitted 
(except for the time, in seconds) and may be assumed to be any consistent unit system. The 
main mass, mi = 2 is connected with a spring, ki = 2, to a fixed point and n = 150 spring-mass 
pairs, kj, ruj, ; = 2, . . . , n + 1, are attached to it. As in [19], the individual natural frequencies, coj, of 
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Figure 2. Nominal response of the main mass to a unit force in the {a) frequency and [b) time domains, with and without the 
additional attachments. 



the smaller attached systems are chosen randomly and uniformally within a range [0, ^^max = 5]. 
The individual masses are randomly drawn from an exponential distribution with expected total 
additional mass, Mj- Therefore, they have probability density, pimjU/Mj) = exp(— (myn/Mj)) for 
all ; = 2, ...,n + l and the individual stiffness coefficients are kj = mjCoj. The actual values of 
the masses and springs, my and kj can be found on the ;-th line of the respective Electronic 
Supplementary Material Files ESM_mn and ESM_omegan, for ] = !,..., 150. Further viscous 
damping, cy = r]kj, is added to each spring, ; = 1, 2, . . . , n + 1 so that the viscous damping matrix 
is C = 7]K for some damping factor 77 = 5 x 10~^. The chosen dimension of the model is such 
that the ensemble of additional masses is not in any asymptotic form of complexity in the sense 
that the n = 150 additional masses cannot really be considered as a good approximation of the 
case of a smooth modal density, notwithstanding the fact that one resonance is isolated. While 
this asymptotic form of an higher complexity model would not pose particular problem to the 
proposed framework, it is not required either. 

The interest is in the transfer function, ^1(0;) = xi{o)) = eJx{oj), of the main mass, mi, owing to 
a unit force, f = ei = [1 0 ... 0]^, as well as in the impulse time response, g^{t) = x^it), of the same 
mass. The evaluated corresponding nominal, 'non-averaged', functions are presented in figure 2. 
Some of the difficulties arising when the complexity of the system increases are illustrated on this 
benchmark problem and now highlighted. First, it is clear that the effect of the additional masses 
cannot be neglected or simply integrated in the properties of the main system. In the frequency 
domain, two main parts of the response are identifiable: a single well-isolated resonance in the 
6 <(jo <7 range and an apparent blurred resonance pattern made of a combination of several 
resonances within the range 0 < a; < 5. While the latter roughly looks like a main individual 
resonance perturbed by a multitude of minor disturbances, such simplified vision is not really 
the case, since it is really the combination of all individual resonances that creates the blurred 
pattern. It is also not clear, de visu and qualitatively, if just one or several equivalent blurred 
resonances are hidden in the pattern. The situation is similarly not that simple in the time 
domain: the impulse response of the full system exhibits some kind of beating phenomena where 
its general magnitude appears to first decrease, i.e. be damped, up to about f = 50[s] to only 
increase back up to f = 150[s]. Similarly, the response at lower times exhibits a high-frequency, 
lower magnitude oscillation, that is superposed to a signal with wavelength equal to about 10[s]. 
Nevertheless, although the impulse response is somewhat complicated, it appears to have some 
hidden simplified pattern, made of only a small number of components. This situation exhibits 
several typical characteristics of the problems one is exposed to when working through various 
frequency ranges. The averaging framework proposed here offers an approach to extract such 
important components. Before discussing this in more detail, the use of two alternative existing 
approaches is presented. 
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Figure 3. Comparison of the magnitude of the nominal (exact) response, (f), to its approximation by modal truncation when 
[a] keeping 1, 2, 5 or 10 modal pairs with the highest modal density magnitude, and [b) keeping 1, 10 and 50 most important 
modal pairs — as defined in [45] — of the attached masses. 



Alternative approaches. An apparently evident alternative strategy to model the response of 
the system would be to consider the modal components of the system that contribute most to the 
response. Looking at the modal magnitudes abs(e]^0y^y^ei) of ejx found from modal expression 
(3.5), as such a criteria, one can truncate the modal representations at a number of pairs of 
conjugate modes (since the response in time is real, each coj comes in pair with an eigenvalue 
with opposite real part, —?il(a)j) + i^{a)j)). The magnitude of the impulse responses are presented 
in figure 3a for an increasing number of 1, 2, 5, 10 modal pairs whose reference eigenvalues, 
sorted by decreasing modal magnitude are coj = {0.5634 + iO.000079, 0.5173 + iO.000067, 6.2552 + 
iO.009782, 0.6022 + iO.000091, 0.4823 + iO.000058, 0.7754 + 0.00015i, 0.6191 + 0.0000961, 0.3755 + 
0.0000351, 0.7165 + 0.000131, 0.8323 + 0.0000241, . . .}. While the predictions obtained through 
modal truncation are roughly similar to the actual impulse response, they are not particularly 
precise. For example, the time and value of the first local maximum, {t,g^{t)) = (0.3800[s], 0.1126), 
of the impulse response is estimated as (0.3350[s], 0.0805), when 10 pairs of modes are used. This 
is a 13 and 40% relative error. Another issue is that it is difficult to assess the quality of these 
approximations if information on the other modes is not used. Other approaches have been 
proposed over the last one or two decades with the objective of reducing model dimensions 
such that the reduced model preserves information of the original system, in particular in the 
time domain. Among the several options that have notably been proposed and explored by 
Barbone et al. [19], it is worth noting a connection of the current work with the smooth (mass) 
modal density function used to model a high-density discrete spectrum. Although starting from 
a different point of view, the present frequency averaging framework can indeed also be seen 
as offering smoothed frequency functions. Further comparison between the two points of view 
might therefore bring additional insight. Further efforts have also been pursued to identify and 
extract the most important modes of complicated subsystems, also with a focus on the quality 
of the time responses of coupled systems. The impulse response predictions obtained through 
the OMR algorithm of [45, Box 1, page 1669]^ are presented in figure 3h for reduced models 
of the attachment of dimensions 1 to 50. While such approach is better targeted than other modal 
truncation methods to the evaluation of impulse responses, its accuracy is still relatively limited 
for the amount of information preserved. Note that the possibly less precise predictions compared 
with those obtained via the previous modal truncation method may be explained by the fact that 
OMR operates without information on the actual system to which the complicated ensemble of 
attachments is connected. The predictions of both sets of impulse responses of figure 3 should 

^Note that the reduced damping matrix was defined as r] times the reduced stiffness matrix. 




be compared to the results presented in figures 4-6. It is seen that an accurate estimation of 
the average can be achieved with only a few equivalent modes and that it can provide extremely 
accurate time domain predictions. A practical implicit approach for evaluating such equivalent 
modes is through the rational Krylov projection method proposed in [31]. 

(b) Evaluation of the frequency moments as computational goals 

Rather than the full, detailed response as presented in figure 2a, it is the three averaged functions 
of order 1 and 2 that are directly targeted in the proposed framework. Such computational 
responses of the system are presented in figure 4 for three different values of the averaging 
width, a. Varying the value of the averaging width allows to look at the response with various 
interest in details, as if using a magnifying glass of various strength. For the largest value, a = l, 
all irregularities of the first blurred resonance are smoothed and it appears from looking at the 
averaged response in (a) that, in some way, two main resonances are present in the response. 
At this value of the averaging or magnifying strength, a = l, the average itself does not allow 
to differentiate between a single, clear resonance and a blurred one. Such information is however 
provided by the second-order averages: the ratio between the square root of the averaged square 
of the magnitude, in (g), and the absolute value of the average, in (a), is about 60 at the first 
resonance, while it is only about 8.5 at the second resonance. This information is sufficient to 
identify the region with higher modal density, and if desired, the averaging width can be refined. 
Doing so reveals more details of the response as seen in (b) where about 10 distinct resonances 
are apparent. Pushing further the identification of details by decreasing the averaging width to 



a = 0.01 leads to a much less regular response in (c). This lightest touch average is however notably 
smoother than the full detailed solution of figure 2: the minor resonances that are not important 
in the response have been smoothed out. 

The averages thus provide qualitative information about the location and density of 
resonances. As such, they are useful by themselves and one can use them to refine the analysis, 
i.e. to decrease the averaging width, in regions of interest. The effect of averaging can also be 
characterized in the time domain, as now described in §5. 



5. Time domain analysis and impulse response of the averaged responses 

In this section, the interpretation and effect of the frequency averaging in the time domain is 
discussed. Constant and variable averaging width, a{a)), are studied in turn. 



(a) Impulse response from frequency averaging with constant width 



Properties of convolutions can be used to evaluate time responses when the averaging width is 
constant. This is described and illustrated in the two following sections, in particular in the case 
of Gaussian averaging. 

(1) Description of the approach to evaluate time responses 

First, for a constant a{a)) = a, averaging the full response, x{co), is equivalent to evaluating its 
convolution with the averaging window or probability density function j)a{/^co). The inverse 
Fourier transform of the average x{o)) is therefore equal to the product of the exact impulse 
response by the inverse Fourier transform of paiAco). The latter is called the characteristic function 
of the probability density function, pa{-)f in statistics. 

Specifically, considering the Fourier transform, x(f ), of the frequency average transfer function 
vector, x{(jo), one has in general the following filtered impulse response: 
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(5.1) 



Since a{a)) = a is constant, the change of variable co^ = (o -\-u gives the product of the impulse 
response of the system by the characteristic function of Pa{-)/ i-e. 
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(5.2) 



(5.3) 



Equation (5.3) offers a strategy to evaluate the time response of the system, which does not 
require the knowledge of the details of the full system, but only its average response. The 
proposed strategy, valid for times at which the values of the scaling function J^^ pa{u)e~'^^^ du 
are large enough, is illustrated in figure 5 for the real Gaussian case. It has three main steps: 



(i) evaluate the frequency average, x{(jo); 

(ii) evaluate its inverse Fourier transform, x{t) and 

(iii) scale this time function, using equation (5.3) in general and (5.7) in the particular real 
Gaussian case, to estimate the impulse response, x{t). 

Each of these operations has the option to be evaluated in various manners. For example, 
while the frequency average could be evaluated by using the modal decomposition used in the 
derivation in this paper, an efficient matrix alternative based on rational approximations and 
Krylov subspaces that has been proposed by Lecomte [31] allows by-passing any modal analysis 
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Figure 5. Comparison of the exact impulse response with its 'Averaged' estimation through frequency average, inverse Fourier 
transform and scaling as described in {g). Three averaging width, a = 1, 0.1 and 0.01 rad Hz, are considered in plots {a), {b) and 
(c), respectively, and the corresponding relative errors are presented in plots (cf), {e) and [f). The inverse Fourier transforms of 
the discretized transfer functions were evaluated through an explicit trapezoidal integral. The real Gaussian scaling functions 
affecting the impulse response when the inverse of the frequency average is taken are presented in (/?). 



and has been successfully applied by the author to large systems in this context of evaluating 
time responses as proposed here. Similarly, the inverse Fourier transform could be evaluated by 
considering individual modes of a modal decomposition or numerically by sampling the average 
in the frequency domain. The evaluation through modal decomposition can also be applied to the 
precise approximate reduced models of [31]. 

Real Gaussian case. In the particular real Gaussian case when p^iu) = Pa (u), the Fourier 
transform of a real Gaussian averaging function is another real Gaussian function: starting from 

^-u^/{2a^)^-iut^^ (5.4) 



iTva^ 



one can complete the square in the exponential argument and use Cauchy integral 
theorem to find 



^i^\i/)e-^"^ du = 



1 ^-{u/{V2a)+ita/V2f^-fay2^^ (53) 



Ina^ 
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The Gaussian filtered impulse response is thus 

x(g)(f) = x(0e-^'"'/2. (5.7) 

This shows that Gaussian frequency averaging for constant a{co) = a preserves the nominal 
impulse response at time zero, i.e. x(0) = x(0) and that the smaller the value of a, the better is 
x^^\t) an approximation of x(t) at small times. 

Note finally that the current discussion provides a description of the eigenfunctions of the 
Gaussian averaging operator. This is detailed in appendix B in the electronic supplementary 
material. 

(ii) Illustration of the approach to evaluate time responses 

The estimations of figure 5 were evaluated through a modal decomposition and the discretization 
of the frequency average in a wide but finite range of frequencies. The relative error on the 
estimated impulse responses is about 10~^ or less in most of the time ranges corresponding to 
each of the averaging widths. The increase in the error at the end of the two first time ranges, i.e. 



the intervals (0,5)[s] for a = l and (0,45)[s] for a = 0.1, is due to two reasons: first, the precision 
of the evaluation of the average and, mostly, of the inverse Fourier transform, and, second, the 
limits of the finite precision arithmetics. Both the evaluation of the average and inverse Fourier 
transform can be made much more precise, particularly by using the method of Lecomte [31]. 
The finite precision arithmetic remark corresponds to a more fundamental fact: since the impulse 
response x(t) has to be scaled by the inverse Gaussian e^^^^/^ to estimate x(^), this estimation will 
be marred by round-off effects if x(t) is represented in finite precision. The estimations for smaller 
values of t and the scaling factor remain precise as is shown in figure 5 where the vertical dotted 
lines correspond to values of this factor equal to ^ = 10^. If an engineer, who would for 
example want to design a structure that is only moderately affected by a shock or turbulence, 
is only interested in the maximum and minimum values of the impulse response during a given 
time period after some impact, as in (a), he or she can thus define the maximum averaging width 
possible that provides accurate response in the smallest time range of interest. The discussion 
of this section provides some explanation of the Gaussian frequency averaging process in that 
the smoothed out details and irregularities of the transfer functions correspond to details of the 
response at larger times. Averaging with other distributions, as those discussed in §3e, leads to 
different but similar interpretations since the characteristic functions of these distributions — their 
inverse Fourier transforms — has other properties. 

Nevertheless their being smoothed out, information about the details of the transfer functions 
does not disappear in the averaging process, at least when the second-order averages — averaged 
square and variance — are considered. The information about the total transient energy, including 
at larger times, is indeed preserved and retrievable. This will be discussed in §6. Before that, the 
case of variable averaging width is briefly discussed. 

(b) Time domain analysis using a variable averaging width 

It was mentioned in §2^ that the low-frequency focus corresponds to responses that have zero or 
small variance. The discussion of the previous section provides the further interpretation that 
their deterministic character can be understood as well defined — rather than blurred — transfer 
functions and consequently refined impulse responses at small times. The alternative high- 
frequency focus corresponds instead to statistical, blurred transfer functions and a more and more 
imprecise notion of energy at larger and larger times. It is now stressed how the focus of the 
analysis can be qualitatively pushed into the direction of a more deterministic character, by 
pinpointing frequency regions. 

In general, since the averaging width, a{a)), can be frequency dependent, it can indeed be used 
to gear the analysis into different focus in different frequency ranges and provide a smooth 
transition from low- to high-frequency as shown in [31]. Another application, illustrated in 
figure 6, is presented here. The frequency average response is by itself an approximation of 
the exact response, in the sense that the averaging process may allow to get rid of unimportant 
features of the response. As noted in the previous sections, this notion of unimportance is related 
to an accurate time response of the system at small times and a small frequency variance. Steps 
for reducing the frequency variance may therefore lead to a focus on more important features of 
the system response. 

The analysis presented in figure 4 for a = 0.1 is examined again. The corresponding inverse 
Fourier transform of the average is presented in figure 6e and is seen to rapidly lose precision for 
increasing time. This can be compared to figure 4h, where the Gaussian scaling was applied. Since 
the frequency variance is relatively the largest in the region oxl, reducing the averaging width, 
i.e. refining the response, in this lower frequency region may consequently achieve a generally 
better approximation in the time domain. 

A constant refined value oi a = 0.004 is chosen in the interval 0.2< co <1 (rad Hz) and is 
connected to the initial a = 0.1 everywhere else through a linear variation at the edges co e (0.1, 0.2) 
and (1,1.1) (rad Hz). With this choice, the average reveals a number (about 30) of resonances 
in the refined interval and globally reduces the variance there. As can be seen in plot (/), this 
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Figure 6. Comparison of the average, variance and impulse response functions for a constant, onqml, averaging width a = 
0.1 rad Hz and an averaging width that has been r^f\n^6 in the interval (0.1,1.1) rad Hz. The approximate and exact impulse 
responses are, respectively, the inverse Fourier transforms of the exact transfer function and that of its average. Both were 
evaluated through an explicit trapezoidal integral as for the results of figures. 



significantly increases the quality of the approximate response over a larger region of time and, 
while qualitative, this study illustrates again the link between averaging and precise response at 
small time. Attention now turns to the imprecise or energy part of the responses that are assessed 
by the second-order averaged square and variance. 



6. Repartition of the energy in the system response 

It is shown in this section that the total transient energy, which transits through part of a 
system during impulse, can be retrieved from a combination of its first and second frequency 
moments obtained with constant averaging width, a{a)) = a. Equivalence theorems of two kinds 
are first highlighted. Exact expressions of the energy in terms of the frequency moments are then 
presented and discussed. It is notably shown how the energy is partitioned in terms of the first 
and the second moment contributions. 



(a) Energy equivalence theorems 

The evaluation of the total impulse energy through averaged functions is based on theorems 
of equivalence of the energy in the three domains, namely the domains of time, frequency and 
frequency average. The known equivalence between the time and frequency domains is reminded 
in the next section while the equivalence between the nominal and average frequency domains is 
introduced in %6a{ii). 

(i) Time and frequency energy density equivalence 

The first kind of equivalence between the frequency and time densities is provided by Plancherel's 
theorem below, whose proof is provided in appendix A (in the electronic supplementary material) 
for completeness. 



Theorem 6.1. Using the notations of equation (2.2) 



X = 



x{t)x{tj^ dt = 



In 



x{oo)x{a)f^ do). 



(6.1) 



Its application to time derivatives gives the expressions needed for the evaluation of the kinetic 
energy, e.g. as follows. 



Corollary 6.2. Using the notations of equation (2.2) 
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(6.2) 



(ii) Frequency and frequency average energy density equivalence 

The second kind of equivalence concerns the frequency integrals of functions of the responses 
and of their frequency average. In particular, the following theorem for the products of responses 
is used for the evaluation of the total potential energy in §6^. 

Theorem 6.3. The integrated quadratic term defined in equation (6.1) is equal to the integrated 
frequency averaged quadratic term, i.e. 
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x{a))x{a))^ da)= — ^_x{a) + /^a))x{a) + Aco)"] dco = X, (6.3) 

■) J— oo 



where E[x(a; + Aa;)x(a; + Ao;)^] denotes the frequency average on the frequency shift Ao; for any 
probability density function, paiAco), and constant standard deviation, a, 



E[x{cD + Aa;)x(a; + Ao;)^] = x{co + Aa;)x(a; + /^cof^pai/^co) d/^co. 



(6.4) 



Proof. Substituting the expression of the frequency average, considering the Constance of a, and 
changing the variable co' = co -\- /^co, one finds X^J^^ {{l/2n) ^'^^x{co')x{cD')^d(j)'\pa{/^co)d/^(jL> 
and one can then integrate successively with respect to co' and dAo; to find X. ■ 

The integrals match similarly when the time derivatives of equation (6.2) are considered, i.e. 

Theorem 6.4. The integrated term defined in equation (6.2) is equal to its integrated frequency 
average, i.e. 
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E[{co + Aa;)^x(a; + Aa;)x(a; + Aco)"] dco = D (6.5) 



for constant averaging width a{o)) = a. 



The proof is similar to that of theorem 6.3 and this property is used for the evaluation of the 
total kinetic energy in §6c. 

These and similar equivalence properties can be applied to evaluate the energy terms of a 
system, just from the knowledge of its frequency average first and second moments targeted in 
the proposed framework. It is now demonstrated how this can be done. 



(b) Repartition and retrieval of potential energy 

The particular case of a typical dynamic stiffness matrix A{o)) = (K + icoC — op-M) is considered 
for illustrative purpose. Examining a transient response of the system, since the potential energy 
of the whole system at any time t is U{t) = (l/2)x(t)^Kx(t), its integral over time is equal to LJ = 
U(t) dt = (l/in) xHHRxM d^o. 



Full system energy. Using theorem 6.3, and denoting the element of K in /cth row and ;th 
column, one has that 
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E[xj{co + Aco)xk{co + Ao;)*] do; 



(6.6) 
(6.7) 



This expression can further be expanded in terms of the average and variance of the response by 
using equation (2.14) which gives 
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or, in compact matrix form, 
1 
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tr 



K 



var[x((X))] d(X» 



x(a;)^Kx(a;) do; 



(6.8) 



(6.9) 



where tr [.] denotes the trace of a matrix, i.e. the sum of its diagonal coefficients. 

Component energy. The integrated potential energy of a component of the system can similarly 
be expressed in terms of frequency averages. For example, if the system is decomposed into two 
components, distinguished by indices 'A' and 'B', such that 



x(a;) = [VA Vb] 
then the energy of the component ; = A or B is equal to 



Xa((^) 


and 


xa(<^) 


_Xb((X))_ 




_xb(<w)_ 



= [Wa WbI'^xH, 



(6.10) 



(i, = ^ jtrj^K^;,- ^ var[xy(c<;)] dw^ + xy(c<;)"%xyHdc<;|, (6.11) 

where var [xy(a;)] = wTvar[x(a;)]W*, Xj{co) = \N^ x{co) and K/y = V|^KVy. Similarly the integrated 
potential energy of the coupling between the two components is equal to Uab + ^ba where 
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tr 



cov[x]^{o)), Xj{o))\ do) 



Xj(oj)^KjkXk(co) doj 



(6.12) 



for Kjk = yfKYk and j,k = A or B. 

For example, the integrated potential energy of a component made of a single degree of 
freedom, say Va = ey where ey is the unit vector with zero components except for 1 at its yth 
component, is 



(var[xy((X»)] + \Xj{(jo)\^) dco 
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■ 471 



E[\Xj{co-\- Aco)\^]dco. 



(6.13) 



Real Gaussian case. In the real Gaussian case, the values of the integral function S^^\co, a, coj) were 
given in equations (3.11)-(3.13) and those of Q^^\(i>,a,coj,co]^) can be found in equations (3.17)- 
(3.25), with appropriate substitution of the right expression of the S^§^ function, based on the sign 
of the imaginary part of its last argument (i.e. the sign of ^{ojj) and ^(oj^)). As discussed in §8, 

the average vector or scalar functions, yP^{(jo) and xP{co), can also be directly evaluated without 
requiring any modal analysis and therefore, by bypassing any evaluation or use of S^^\co,a, coj). 



(c) Repartition and retrieval of l(inetic energy 

The integrated kinetic energy can also be expressed in either time or frequency domains. The 
expressions of the variance and products cannot be used directly owing to the presence of the clP- 
term in the right-hand side of the last equation. The integrand must first be modified. Starting 
from equation (3.5), one finds 



i=l,...,N ^ J ' ;=1,...,N 



1 



(6.14) 



which allows to evaluate the corresponding frequency average in terms of known integral 
functions S{o}, a, coj) and Q{co, a, coj, co]^). Indeed, for a general pdf pa{^(^) 

E[{a) + Aa)fx{a) + Aa))x{a) + Aco)^] 

(oj + Aco)^x{oj + Aoj)x(co + Aco)^pa{Aoj) dAco 



k^l,...,N j=l,...,N 

X [1 — (jL>jS{(jL>,a{co),(jL>j) — co^S{co,a{(jL>),co^) + cjL>jCo^Q{co,a{co),coj,cjL>]^)]}. 



(6.15) 



This modal expression of the second moment can also be expressed in matrix form, by using 
expression (3.10) of the averaging matrix H and the general discussion of §3a. Basic algebra, 
together with the substitutions #(&>) = H^^\co)i^^\ veiY[x^^\co)] = E[x^^\co + Aco)x^^\co + Aco)^] - 
#(a;)#(a;)", F^{Af)-H^^^ = 0 and F^(Af)-^A^Q^ = [0 - 11]^, gives indeed 

E[{co + Acofx{co + Aco)x{co + Aco)^] = d{co)d{oj)^ + var[d(a;)], (6.16) 

where d{a)) = ia;x(a;), d{a)) is its frequency average, retrievable from the second part of the vector 
#(a;)andvar[d(a;)] = [0 - l]var[x(l)(a;)][0 - if. 

In the typical particular case of the previous section, i.e. A{co) = (K + icoC — o/M), the integrals 
of the kinetic energy, T = J^_^(l/2)(9x(f)/9f^)M(9x(f)/9f) dt, can thus also be expressed in terms 
of computable averaged frequency densities thanks to these expressions. 

Full system energy. Through theorem 6.3, one knows that the kinetic energy of the whole system 

is 



k=l,...,N ;=1,...,N ^ 
k=l,...,N ;=1,...,N ^ 



oj^Xj(oj)x]^(oj)'^ doj 

) 

E[{co + Acofxjico + Aco)Xk{co + Acof] dco 



(6.17) 
(6.18) 



This can further be expanded in terms of the average and variance of the response, by using 
equation (6.16), and it can finally be written in a compact matrix form similar to that of equation 
(6.9), i.e. 



T = T-- 



var[d((X»)] dco 



d(a;)"Md(a;) doj 



(6.19) 



Component energy. The kinetic energy of two individual components distinguished by indices 
'A' and 'B', generally such that 



d(a;) = [VA Vb] 



dA(<^) 


and 


dA{oj) 


_dB{o))_ 




_d^{oj)_ 



is then 



tr 



var[dy(tti)] dco 



+ 



:[Wa WefdH 



dy(ft))"M^ydy(ft)) dft) 



(6.20) 



(6.21) 



where var[dy(a;)] = [Wy]Tvar[d(a;)][Wy]*, dj{a)) = [Wj]^d{a)) and M^y = Vj^MVy, and ; = A or B. 
Similarly, the integrated kinetic energy of the coupling between the two components is equal 
to Tab + Tba where 



cov[dk{co), dj{a))] dco 



dj(oj)^Mjkdk(oj) dco 



(6.22) 



for;,/c = A or B. 

The kinetic energy of a component made of a single degree of freedom, say the ;th one such 
that Va = ey, is thus the following sum of two terms: 



Ta- 



1 

4^ 



veiY[dj{oj)] dco + |dy(a;)p dco 
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E\dj(co)\^ dco (6.23) 



where all the averages of the scalar response dj{co) may be evaluated as discussed. 

Real Gaussian case. As for the potential energy discussed in §6^, in the case of a real Gaussian 
average, the modal expressions only necessitate the functions S^^\co,a,coj) and Q^^\co,a,coj,co]^) 
and, as discussed in §8, the matrix expression of the average can be evaluated without requiring 
any modal analysis. 



(d) Retrieval of input, dissipated and other energy terms 

Other energy terms can be similarly evaluated exactly by integrating frequency averaged 
quantities. For example, the energy dissipated by viscous damping necessitates the evaluation 
of a term {l/ln) J^^ cox{co)x{co)^ dco = {l/ln) J^^ E{{co + /^co)x{co + Aa;)x(a; + Aco)^] dco whose 
second expression can be evaluated using equation (6.14). 

The integral of the transient input power can also be evaluated in terms of the averaged 
expressions. For example, in the case of an impulse force, its frequency density is constant so 
that the total input energy is 

{i^{t))^x{t)dt-- ^ 



£(in) 



In 



f^x(a;) dco 



and this term can again be expressed in terms of a frequency average 



£(in) ^ £(in) ^ _fT 

In 



E[x(a; + Ao;)] do; = — 

. In 



x{co) dco. 



(6.24) 



(6.25) 



(e) Averaged square and averaged product of responses 

The averaged square and product of responses did not appear in the expressions of the energy 
terms in the previous sections. They are however an important feature of the response at higher 
(mid or high) frequencies and are therefore an integral part of the computational goals of the 
proposed framework. Without these notions, a pure energy analysis indeed loses information on 
the phase of the transfer functions that may possibly be critical for some applications. Although 
these terms should not be overlooked, they may of course be omitted in situations where they 
are negligible, or in asymptotic forms, in the same way as the variance may be omitted in the 
asymptotic case of a low-frequency, deterministic, analysis. 



7. Time domain analysis of the covariance and averaged products 
of responses 

In this section, an interpretation of the second frequency moments is provided in the time domain. 
The discussion is similar to that of %ba where it was shown that, since frequency averaging with 
constant width is the convolution of an averaging function by the exact transfer function, its 
inverse transfer function is the product of the exact impulse response by the inverse Fourier 
transform oipa{^<j^)- Similar properties exist here with regards to double or two-dimensional inverse 



Fourier transforms of the second-order averaging terms, that is of the covariance and averaged 
products of responses. These properties are now detailed, with a particular focus on the real 
Gaussian case. 

Starting from the definition of the inverse Fourier transform of equation (2.2), one finds 



CSq[XA((^A),XB((X)B)]^ 



xa(<^a + u)xb{ojb + uYpa{u) du 



xa(^a) 



XBaB)^csq[e-i"A^A^ e-io^B^B] d^A, (7.1) 



where the integral in csq[e"^^A^A, e"^^B^B] = e^^^^^^ e^^^^^^ pa{u) e-'""^^"^^^^^ du needs to be 
evaluated for particular functions, pa{-)- Following the same derivation, for the case of the power 
term, one finds 

E[xA(a;A + Ao;), x^icoB + Aco)] = xa(<^a + w)xB(a;B + uf^paiu) du 



xa(^a) 



XB(fB)csq[e-i"A^\ e+i"^^^] d^B d^A- (7.2) 



Real Gaussian case. In the real Gaussian case, pa(.) = Pa K^r orie can again use equation (5.6), 
which gives 

p^^\u)e-'<'^+'^^du = e-^'^+'^'>"'"'^ (7.3) 

and therefore 

CSq^S)[xA(a;A),XB(tf^B)] 

^oo / poo J- . \ 

XA(f A)xB(fB)^e-[('*+^'') " /^l e-'""''" dfB e-'-^^^* dfA (7.4) 

J—OO \J— OO ^ / 

and 

E[XA(a;A + Ato), XB(a;B + ^0))^\ (7.5) 

-OO / poo , . \ 

XA(f A)xB(-fB)^e-[('A+''') " e-'""'''' dfB e-''"*'* dtA, (7.6) 

J—OO \J— OO ^ J 

where, in the last expression, the argument of x^ is negative. The variance derives directly from 
this expression. 

Each of the two expressions (7.4) and (7.6) is a two-dimensional Fourier transform in the 
plane (^A/^b) of the product of an outer product of Xa and Xg by the weighting function, 
^(^A/ ^b) = e~[^^^+^^^ ^ While similar expressions would exist for other averaging distributions, 
the weighting function in the real Gaussian case, is remarkably also a Gaussian. Furthermore, 
since both impulse responses Xa and Xb ^^e zero for negative value of their arguments, in the case 
of a causal system, and since Xg has a negative argument in the expression of the second average, 
these inverse Fourier transforms of the average functions result in functions that are each non- 
zero only in a single quadrant. In both cases, the outer product of the impulse responses is scaled 
by the Gaussian ridge function, r(f A/ ^b)/ that becomes tighter and tighter for increasing values of a. 
The maps of the non-zero components of these functions are illustrated in figure 7 for the transfer 
function gi, which was introduced in §4. 

Different behaviours are visible for different averaging widths. For the smallest, a = 0.01, the 
averaged products and covariance are very close to the non-averaged outer products of impulse 
responses. Simplifications however occurred since irrelevant details have been smoothed out 
of the response. This corresponds to a low-frequency behaviour. At the largest, a = l, the cross- 
influences are negligible while the energy information is preserved. This is the high-frequency 
behaviour. For the intermediate value, a = 0.1, additional information about the coupling between 
responses at different times, i.e. phase in the frequency domain, is partially present. Observing 
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Figure 7. Illustration of the two-dimensional inverse Fourier transform of the second frequency moments as described in (g). 
Three averaging width, C7 = 1, 0.1 and 0.01 rad Hz, are considered as well as two kinds of second frequency moment terms: 
plots {a), ib) and (c) are the maps of non-zero values of the inverse Fourier transforms of the frequency averaged products, 
csq^9^[^i(a;A)A^i(<^B)]A whereas plots (cf), (e) and (0 are those of the power term E[^i(a;A + A(jo),g]{(jo^ + A(jo)\ both 
for a real Gaussian average. They are respectively labelled 1F(E[^i^|])' and 1F(E[^i^{^])'. While the specific transfer function 
of §4 is used here for illustration, the same Gaussian ridge scaling functions, r(ti\, t^) presented in (/?), would affect the 
product of the exact impulse responses of any general system similarly averaged. 



both the covariance and averaged products provides a rational and quantitative tool with 
which to assess the actual behaviour of the dynamic system. As such, the quality of the existing 
methods to deal with the frequency methods can be analysed and compared by evaluating how 
well they preserve the actual frequency average, covariance and averaged product characteristics 
of the actual system. New methods can also be developed based on the actual important system 
features identified by the averaging process. 



8. Efficient evaluation 

Besides the asymptotic approximations, which provide evaluation of the averages that are 
accurate and economic in particular situations, the question of evaluating the averages efficiently 
and precisely in a general context is a relevant concern. As mentioned in §3e, a solution has 
been offered in [31] for the first moment averages, through the precise expansion of distribution 
or weighting functions, pa{AQ), into fast converging rational expressions. Consequently, the 
exact averages are approximated in terms of fast converging rational functions of the poles of 
the system. These rational functions may be expressed in matrix form and evaluated exactly 
through stable Krylov projection methods, using only solutions of a system with added damping. 
Combining the projection subspaces at only a few interpolation frequencies can then provide a 
precise evaluation of the average function over a whole frequency range. This approach has been 
shown to be extremely efficient and able to provide the Gaussian averaged response in a whole 
frequency range, and with frequency varying averaging width, with only a few responses — as 
little as two — of the full system. This has particularly proved valuable to obtain the response of 
the system through a single reduced model of the system that gives different low- to high-frequency 
focus through different frequency ranges. The same strategy of mixing rational expressions, 
interpolation, and model reduction is applicable to other weighting functions that may or may 
not be originally expressed in rational form. In terms of asymptotic approximations, different 
approaches can be used, for example with rational expansions being used at zero, intermediate 
or infinite frequency. 



The question of the efficient approximation or representation of the covariance and averaged 
product is of essential importance for the fundamental understanding of the mid-frequency region. 
The framework and analysis of this paper offers a new point of view, notably in that the pairs 
of plots {a)-{d), (h)-{e), {c)-{f) in figure 7 expose that the relevant patterns are the product of an 
outer product of functions in the (^A/^b) plane and another outer product of a Gaussian with a 
constant functions in the (f a + ^B/ ^a — ^b) plane. The consideration of the averaged product of 
responses additionally to the covariance is important in assuring this representation, as well as 
in preserving information on the phase, as previously mentioned. It is evident from the analysis 
and plots that it can, in theory, be retrieved from the covariance or energy information evaluated 
with a single averaging width over the whole frequency range. 

The question of optimality of the approximations of the responses as in [19] can be reframed 
and possibly extended in the proposed averaging framework. 

9. Conclusion 

A solution framework for dynamic systems has been proposed that targets directly the frequency 
average of their responses and their higher, second frequency moments. This approach is 
applicable exactly, independently of the system's modal density or modal sensitivity so that 
a smooth transition from one region of the frequency spectrum to the other can be achieved 
by merely tuning an averaging width function, a{(jo). This parameter can be interpreted as the 
typical standard deviation or the characteristic bandwidth of a distribution or weighting function 
Pa{Aa)). It describes the desired level of frequency precision or uncertainty of the solution at any 
frequency. Asymptotically zero values of the parameter correspond to a deterministic or low- 
frequency solution with zero variance, while larger values correspond to high-frequency solutions 
that can be understood in a more energetic or statistical sense. The transitional mid-frequency 
region is smoothly supported and no distinction or boundaries between the regions is necessary. 
The frequency averaged power or energy of the system, which are statistics that are usually 
considered at high-frequency, are found as function of the frequency average and of the averaged 
square of the absolute value of the response. It has been explicitly stressed that this fact considered 
in [38] for the case of a uniform distribution can be applied for a general weighting function. 
The general explicit analytical expressions are provided and demonstrated here for the case of a 
real Gaussian averaging function. Again, the flexibility in the choice of the frequency-dependent 
parameter, a{co), allows to extract deterministic, statistical or intermediate frequency averaged 
power, with values of the parameter that are respectively small, large or intermediate. Contrary 
to the case of many statistical energy analysis and high-frequency methods, the power or power 
average is not the only result of the analysis in the proposed framework. The consideration of 
both frequency average and variance, and of the tunable uncertainty parameter, a{oj), in a single 
analysis gives a smooth transition from low- to mid- and high-frequency ranges and provides a 
framework in which the usually distinct deterministic and statistical approaches can be seen just as 
particular or asymptotic cases. 

Many existing low-, mid- and high-frequency approaches can actually be placed and examined 
in the proposed framework. For example, deterministic modes, models of low-frequency physical 
components, such as springs and masses, SEA systems, analytical waves or hybrid models can 
be integrated into a single model. The coupling between any of these components within the 
framework must be characterized in such a way that the average, averaged square and variance of 
the responses of the global system are accurately represented. This may be seen as a generalization 
of the fact that the coupling of high-frequency or SEA components must be such that the transfer 
of energy in a whole SEA system is accurately represented and that the coupling between low- 
frequency or deterministic systems must be such that the deterministic transfer functions, i.e. 
their frequency average for zero value of the a{.) parameter, of the coupled system are properly 
predicted. Coupling conditions assuring matching of the averages and variance also offer a 
point of view for the coupling conditions between the deterministic and SEA components of 
hybrid models [46,47]: for the same averaging width, a{(jo), a sub-system might behave in its 



low-frequency regime while another behaves in its high-frequency regime in which case, this 
would be reflected in the average and variance of the response of the components. Low-, mid- 
and hybrid mid-frequency coupling conditions are also particular cases of the more general 
framework coupling conditions. They are perfectly valid within the framework but only under 
certain conditions and for specific ranges of values of the parameter a{.). This is exactly in the same 
manner as the validity of a simplified model of a spring or a mass starts breaking down when the 
frequency increases and the wavelengths become closer to the dimensions of the actual physical 
spring or mass. Other topics of the literature that have not yet been covered here but would 
be worth studying in future works include the averaging in the location of force, measurement 
and interface areas, the treatment of higher moment such as the variance of the general energy 
terms, and the conjunction of frequency and statistical averages. On the other hand, material more 
rarely covered in the literature, that has been studied and highlighted here includes the treatment 
of the averaged product of responses and the covariance terms and the interpretations of the 
averages in the time domain. An efficient approach based on the frequency averages to estimate 
time responses has notably been proposed. 

Besides the fact that many approximate or asymptotic methods can be integrated into a single 
analysis, the framework also provides the advantage of offering a natural environment in which 
they can be better understood and expanded. The precision of various methods to approximate 
low-, mid- or high-frequency components can also be studied in reference to how they affect the 
precision of the average and variance of transfer functions. 
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